Nonequilibrium dynamical mean field simulation of inhomogeneous systems 



m 
o 



Martin Eckstein 1 and Philipp Werner 2 

1 Max Planck Research Department for Structural Dynamics, 
University of Hamburg-CFEL, Hamburg, Germany 
2 Department of Physics, University of Fribourg, 1700 Fribourg, Switzerland 
(Dated: March 20, 2013) 

We extend the nonequilibrium dynamical mean field (DMFT) formalism to inhomogeneous sys- 
tems by adapting the "real-space" DMFT method to Keldysh Green's functions. Solving the coupled 
impurity problems using strong-coupling perturbation theory, we apply the formalism to homoge- 
neous and inhomogeneous layered systems with strong local interactions and up to 39 layers. We 
study the diffusion of doublons and holes created by photo- excitation in a Mott insulating system, 
the time-dependent build-up of the polarization and the current induced by a linear voltage bias 
across a multi-layer structure, and the photo-induced current in a Mott insulator under bias. 
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I. INTRODUCTION 

The study of nonequilibrium phenomena in correlated 
lattice systems has become an active research field due 
to experimental progress on several fronts. In cold atom 
systems, the interaction and bandwidth can be con- 
trolled via Feshbach resonances and the depth of the lat- 
tice potential, respectively, while the effect of external 
fields can be mimicked by shaking or tilting the optical 
lattice!" This allows to investigate quench-dynamics or 
field-driven effects in systems which may be viewed as 
ideal realizations of the simple model Hamiltonians typ- 
ically considered in theoretical studies. On the other 
hand, advances in ultra-fast laser science have made it 
possible to perturb a correlated material with a strong 
pulse and track the time-evolution of the system with 
the (femto-second) time resolution needed to observe in- 
trinsically electronic processesi&£ Such experiments can 
provide new insights into the nature of correlated states 
of matter and may even lead to the discovery of 'hidden 
phases', i.e. long-lived transient states that cannot be 
accessed via a thermal pathway. 

Stimulated by these developments, a growing theoret- 
ical effort is aimed at describing and understanding the 
nonequilibrium properties of correlated lattice systems. 
Given the complexity of the task, much of this work 
has focused on the simplest relevant model, the one- 
band Hubbard model, which describes electrons that can 
hop between nearest-neighbor sites of some lattice with 
hopping amplitude t, and interact on-site with a repul- 
sion energy U . A method which is well-suited to cap- 
ture strong local correlation effects is dynamical mean- 
field theory (DMFT),->£ and this formalism can be ex- 
tended to nonequilibrium systems in a rather straight- 
forward manneri 10 ' 11 Over the last few years, nonequilib- 
rium DMFT has been used in a large number of theoret- 
ical studies of the nonequilibrium dynamics in homoge- 
neous bulk systems, including interaction-quenches j 12 ' 13 
dc-field driven Bloch oscillation o 11 ' 14 or insulator-to- 
metal transition s 15 ' 16 (and the related phenomenon of 
dimensional reduction^) , photo-doping^ ac-field in- 



duced band-flipping, 19 and nonequilibrium phase transi- 
tions from antiferromagnetic to paramagnetic states . 20 ' 21 
While connecting these results to actual experiments is 
difficult because of the idealized set-up in the model cal- 
culations, they have provided important insights into the 
relaxation dynamics of purely electronic systems, and the 
associated time-scales and trapping phenomena. 

One step towards more realistic model calculations 
is to switch from infinitely extended, homogenous sys- 
tems to a description which allows for a spatial varia- 
tion in the model parameters. In equilibrium, the in- 
homogeneous DMFT approac h 22 ' 23 allows, for example, 
to describe some effect of the trapping potential in cold- 
atom experiments ; 24 ' 25 or correlation effects in artificially 
designed heterostructuresj 22 ' 23 ' 26 In a direct generaliza- 
tion of this real-space DMFT to nonequilibrium, one 
would have to store and manipulate Green's functions 
Gij(t,t') which depend on two space arguments i,j and 
two time arguments. Decoupling of space and time is 
no longer possible, neither by introducing momentum- 
dependent Green's functions Gk(t,t') (as in homoge- 
neous nonequilibrium DMFT), nor by using frequency- 
dependent Green's functions Gij(uj) (as in inhomoge- 
neous equilibrium DMFT). The fully inhomogeneous set- 
up would thus require a prohibitively large amount of 
memory for most applications. However, the problem 
turns out to be numerically tractable for a simpler lay- 
ered geometry, which is still relevant for many applica- 
tions. Here one considers a system in which the proper- 
ties can change as a function of the lattice position in one 
direction, while being homogeneous in the d — 1 other di- 
mensions. For example, such an extension allows to deal 
with surface phenomena in condensed matter systems, 
such as the propagation of excitations from the surface of 
a sample into the bulk (which has been looked at recently 
within a time-dependent Gutzwiller approach^!) . In this 
context it is important to mention that pump-probe ex- 
periments often excite only a thin surface layer, such that 
interesting phenomena must be inferred by subtracting 
the bulk-signal, based on some assumptions about the 
penetration depth of the pump pulse. The layer descrip- 
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FIG. 1: Illustration of the layer set-up with N correlated 
layers (full dots), intra-layer hopping t, inter-layer hopping t x , 
interaction U, and chemical potential [i (all these parameters 
can be layer- and time-dependent). The boundary condition 
is given either by coupling to some noninteracting equilibrium 
bath, by a vacuum (no hopping into the boundary layers), or 
by repeating the hybridization functions of the surface layers. 



tion also naturally lends itself to the study of interfaces 
and heterostructuresi 28 i 29 The latter are at present the 
subject of extensive research, and experimental results 
on ultra-fast photo-induced metal-insulator transitions in 
heterostructures have recently been published.— 

In this paper, we discuss and test an implementation 
of the nonequilibrium DMFT formalism for inhomoge- 
neous, layered structures. This formalism is an adapta- 
tion of the equilibrium "real-space" DMFT method de- 
veloped by Potthoff, Nolting, Freericks and othersi 22 i 23 
We discuss the formalism and the techniques used for 
solving the DMFT equations in Sec. lU and illustrate 
the versatility of the approach in Sec. Mil with several 
test calculations involving electric field pulse excitations 
of correlated layers or heterostructures. Section ITVl gives 
a brief conclusion and outlook. 



II. MODEL AND METHOD 

The approximate DMFT treatment of layered struc- 
tures assumes a local self-energy for each layer and maps 
the system to an effectively one-dimensional model sub- 
ject to a self-consistency condition. An efficient strat- 
egy for solving the DMFT equations, which involves a 
partial Fourier transformation of the Green's functions 
with respect to the transverse space directions, has been 
proposed by Potthoff and Nolting^ and a detailed de- 
scription of the equilibrium implementation of this so- 
called "real-space DMFT" can be found in work by 
Freericksi 23 ' 31 Here, we adopt this technique to nonequi- 
librium systems in order to describe pulse-excitations of 
surfaces or heterostructures, as well as transport through 
correlated thin films, using the nonequilibrium DMFT 
formalism! 10 ! 11 



We consider a Hubbard model with N layers, con- 
nected by an interlayer hopping t 1 - , and either "vacuum" , 
"lead" or "bulk" boundary conditions applied to the left 
(n = 1) and right (n = N) surface layers (see Fig. [TJ. 
Here, "vacuum" means no hopping to the boundary layer, 
"lead" means we impose some equilibrium DMFT solu- 
tion in the boundary layer, and "bulk" means that the 
solution on the surface layer is repeated periodically. The 
corresponding Hamiltonian is given by 

N 

H = [ _ C l,n,tr c 3,n,a + £loc,n c i,n,cr C i,n,<7 

n— 1 ija ia 

N 

n— 1 i 
N-l 

+ J2 J2 ( - ^An^n+l.a + h.C.) + b.t, (1) 
n—1 i(T 

where cj creates an electron on lattice site i in layer n, 
U n is the on-site Coulomb interaction in layer n, ei oc ,n is a 
layer- dependent on-site potential, and t" and t 1 - denote 
the hopping within the layers and between the layers, 
respectively. The term "b.t." summarizes the boundary 
terms as described above. All parameters can depend 
both on time and on the layer index, which will mostly 
not be shown explicitly in the following. In the actual im- 
plementation, each layer corresponds to a (^-dimensional 
hypercubic lattice with lattice spacing a, and we present 
results for d = 1. We will later switch to the Fourier 
transformation with respect to the intra-layer coordinate, 
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Ck,n,a- The intra-layer hopping 
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Hamiltonian becomes J2ka e ^,k c k a Ck y n,a, with the dis- 



persion e„ t 



V i" 



^ik{rj —Ti) j a 



External electromagnetic fields are included in Eq. (JTJ 
via the Peierls substitution: We consider electric fields 
E = (En,E^) that depend only on the layer coordinate, 
and let E„ and E^ denote the parallel field component in 
layer n and the perpendicular field-component between 
layer n and n+l, respectively. Units for the fields are 
taken as [t]/ea for E^ and [t]/ea ± for E ± , where [t] is 
a unit of energy, a ± is the spacing between layers, and 
— e is the electron charge. In a gauge where also the 
scalar potential <j> n and vector potential A = (An,A^) 
depend on the layer only, we then have E„ = —dtAn and 
E£ = ~d t A^ — ((f> n +i — 0n)j an d the Peierls substitution 
gives 



^loc,n — Qoc,n T*ni 



(2) 
(3) 
(4) 



where quantities with tilde correspond to zero field. Also 
for fields perpendicular to the layer it is often convenient 
to use a gauge with zero scalar potential. 

Nonequilibrium DMFT provides a set of equations 
for the space- and time-dependent Green's functions 
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Gi t n;j,m(t,t) = -*(7c c J,™^(*) c ],m, CT (*'))- Here t and *' 
lie on the L-shaped Keldysh contour C, and 7c is the 

contour-ordering operator. The notation for contour- 
ordered Green's functions and their inverse operators 
is adopted from Ref. [l3j. The functions Gi jTl; j tm (t,t') 
are obtained from the lattice Dyson equation with a 
local but layer-dependent self-energy £„(t, t'), which is 
computed from an effective impurity model (see below; 
for simplicity we omit a possible dependence of local 
quantities on spin). Due to the translational invariance 
within the layers, one can perform a Fourier transfor- 
mation in the transverse directions and introduce the 
momentum-dependent Green's functions Gfe ;n , m (t, i') — 
-i(TcCk yn ^(t)cl m rJ (t')}. The Dyson equation then de- 
couples for each fc, and one has the following matrix ex- 
pression for the N x N matrices (Gk) n ,m = Gk;n,m, 



which involve the Green's functions G^ for the "chain" 
(Eq. (J5J) ) with site n removed. The Green's functions 
G^ satisfy equations analogous to Eq. (TTU]) . such that 
we obtain for the hybridizations A^ and A R 



k,n 



\L _ ,±* 



9k. 



fc,n— 1 



L ^k,n — l n-l * -1 



a ' - A R 



*ti 



(15) 
(16) 



(Gfc 1 ) 



[idt ~l~ M ^loc,m ^kjn ^m)^m,n 



for layers n = 1,...,N. The boundary conditions read 
A fe>n = ("vacuum") or A fci „ = A fe!lcad ("lead") for 
n = 0, N + l. The "bulk" boundary condition is A§ +1 = 
A§, A£ = Af. Once the A£ n _j and A R n+1 for a given 
layer n have been updated, one computes (Gk) n ,n using 
Eq. (|10p , and determines the hybridization function A„ = 
A n [G„] of the impurity model by solving the impurity 
Dyson equation, 



which is equivalent to the Dyson equation for a one- 
dimensional chain with sites m — 1,...,N. The lo- 
cal Green's function on layer n is then computed from 



G n — — y^(G'fc)„ ;n 



1 



id t + /i - eioc 



A. 



(17) 



Gn - w:T,k(Gk)n,n, and hence we only need the di- The solution of the impurity problem (in the present case, 



agonal elements (Gk) n ,n of the momentum-dependent 
Green's function. These can be evaluated using the fol- 
lowing formulas for the inverse of a tri-diagonal matrix, 



/ z — a\ b\ 

b\ z ~ a 2 b 2 



\ 



we use the non-crossing approximation (NCA ) 3 ' 33 as im 
purity solver) yields an updated G n and £„. 

A self-consistent solution on all layers can hence be 
obtained by the "zipper algorithm" : 23 
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(7) 
(8) 
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t A^ , G n , A „ I A^ , G„ , A„ 



Explicitly, one finds 

(Gk)n,i 



a -1 - A L - A R 



(10) 



Vkli = id t + M - eioc,n - £fe,« - S„, (11) 

where gk,n is the Green's function corresponding to an 
isolated layer, and we have introduced the products 

=*n-lWGL"L-l(*.*')*n-l(*') (12) 

A£ n+ i(t,0 =*n(*)GtUi(*.OC(*'). (14) 



where we start for example with A„ = Abuik, S„ = Sbuik, 
for n = 1, . . . , 7Y, Af = A$ = 0, and then update A£ 
using Eq. (fTS"|) from n — 1 to AT. On the way back, we 
use Eq. (|T6l) to update A^, from n = N to n = 1, and at 
the same time compute G n , A„ and £„ for each of these 
n, and so on. 

Equations ((TU)) - ((T7)) are integral-differential equations 
on the Keldysh contour. Following the strategy outlined 
in Ref. [l8l we can cast these equations in a form that 
can conveniently be handled by numerically stable "time- 
stepping" procedures for the propagation of Green's func- 
tions in real time. Defining the variables = ek,n + 

^k,n-l + ^k,n+l and Z n = [ id t + fi ~ eioc,™ ~ S„] _1 , One 

can write Eqs. pU|) and (fT?) in the form (dropping for 
simplicity the index n everywhere) 



[z- 1 - a] * G k = 1, 

[Z- 1 - A]*G = /, 



I = G h *[Z- 1 -£ k ], (18) 
/ = G*[Z- 1 -A]. (19) 



By summing Eq. (|18[) over k and comparing with 
Eq. JTSJ, one finds G^ = J] fc Cfc * Gfc = A * G and 

Sfc Gfc * Cfc — G * A. We next take the second of 
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Eqs. (IT9|) and multiply from the right with Z. This leads 
to 

[I + G*A]*Z=[I + G (1)J< }*Z = G, (20) 

which we can solve for Z (after having evaluated G^ = 
G * A). Multiplying the first of Eqs. (TH) from the left 
with Z gives 

[I-Z*£ k ]*G h = Z, (21) 

which we can solve for G k . From the first Eq. (flT)|) . we 
also get [I + A*G]*A = Z~ 1 *G*A = J2k Z^ 1 *G k *i k . 
But from the first Eq. (fill), z ^ * G k = I + Cfe * G k , so 

[/ + A*G]*A= [7 + G (1) ]*A = G (2) , (22) 

where G^ 2 ^ = J2 k (£k + £,k * G k * £ fe ). The solution of 
Eq. (HU) yields A. 

To solve Eqs. t|15|) and l|16p. we first compute <jffc i?1 , 
[Eq. (TIT))] by solving 

[J - Z n * efe] * g k , n = Z n . (23) 

With the shorthand notation = (gJT +1 ')„ and 

„ = (G[! l_11 )„, we rewrite Eqs. and as 

[/ - <?fc,n * Cl * * tn-l] * A fe,n = 9k,n, (24) 

[/ - g k , n * ti * ££ n+1 * ti*] * A« n = fffe ,„. (25) 

Equations {U}, 42TJ), (22), ((Ml), (23) and (25), are all 

of the form [J + A] * X = B and have to be solved for X. 
This is an integral equation of the Volterra type, which 
is well behaved and which we solve using the techniques 
described in Ref. [TH. The solution can be obtained by 
successively increasing the maximum time in a step by 
step manner, thereby not modifying an already converged 
solution at earlier times. 

In summary, at a given time-step, we perform the fol- 
lowing calculations in layer n: 

1. For given A„, solve impurity problem (NCA equa- 
tions) to obtain G„. 

2. Evaluate G^ = G n * A„, solve Eq. (20) for Z n . 

3. For each fc-point, 

• solve Eq. (|2"3")l for g k . n , 

• solve equations of the type (j2"5j) and (23) to 
get the new A% n or A£ n , and compute Aj* 
or A£ n from Eqs. (T2) and ([HI) (depending 
on the direction of the sweep), 

• define £ fe ,n = £fe + Afc,«-i + A^„ +1 and solve 
Eq. CD for G fc ,„. 

4. Having obtained and G k%n for all fc-points, cal- 
culate Gn and G^ 2 ^. 



5. Solve Eq. (f22|) to obtain the new A„. 

Then we move to the next layer, where we repeat the 
same cycle, zipping back and forth until convergence is 
reached. Only a few cycles are needed for convergence, 
since a very good starting point is obtained by extrapo- 
lating the Green's functions from earlier times. 

Depending on the application, it may be desirable to 
include a dissipation mechanism which allows to remove 
energy injected into the system by a quench or external 
field. In Ref. HH we have briefly described how one can 
locally couple a phonon bath with given temperature. 
Let us discuss now how such a bath can be incorporated 
into the "zipper algorithm" . In our approximation, the 
electronic self-energy on layer n is the sum of an elec- 
tronic contribution, Ej/[G n ], and of a bath contribution 
Sdiss [G n ] ■ As in the case without bath (Eq. (fl?) ) , S;y [G n ] 
is obtained from the solution of the impurity problem 
with hybridization A n : G n = G„[A„], with 

Gn = id t + n-Vu[G n ]-A n ■ (26) 

The bath contribution is approximated by the lowest or- 
der Holstein-type electron-phonon diagram: 

E diss [G n ] = XG n (t,t')D(t,t'), (27) 

with D(t,f) = -iTr[T c exp(-i f e dtw tfb)b(t)tf(t')]/Z 
the equilibrium boson propagator for boson frequency 
uj and coupling strength A. Therefore, in Eqs. (fT0|) 
and (|llj) . which relate the momentum dependent lat- 
tice Green's function to the self-energy, we have to re- 
place E„ by £[/[G„] + £di SS [G„], or equivalently e k by 

€ k + E(j; ss [G n ]j^ 

In practice, we define Z\ att = [id t + /x — ei oc ,n — ~ 
^diss] -1 (dropping the layer-index n), so that Eq. ([2"o]) 
becomes G = l/O^iatt - A iatt), with Ai at t = A - T, diss . 
We may then repeat the derivation of Eqs. (fIT)) - (j2"2") with 
the substitution A -)• Ai att , Z -t Z\ at t, i.e., given Aiatt 
and G, Z latt is computed from [I + G^] * Z latt = G 

(with GW T = G * Ai a tt), then a new Ai a tt is obtained 
from the solution of [I + G&] * A latt = G {2 \ Finally, 
A = Ai a tt + Ediss is used as input for the impurity solver. 

III. RESULTS 

A. Test of the implementation 

In this work we will consider 1-dimensional layers, and 
use the intra-layer hopping i" = 1 as the unit of energy. 
The equilibrium spectral function for an infinite system 
of such 1-d layers and inter-layer hopping t x = 1 (corre- 
sponding to the usual 2-d Hubbard model) is shown in 
Fig. [2] for inverse temperature (5 — 5 and indicated val- 
ues of U. The impurity problem was solved with NCA on 
the Keldysh contour and the spectra were obtained via 
Fourier transformation of the retarded Green's function. 
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FIG. 2: Equilibrium spectral functions for an infinite system 
of 1-d layers, at /? = 5 and indicated values of the interac- 
tion strength. The DMFT solution has been obtained with 
an NCA impurity solver, which yields reliable results in the 
insulating phase. 



Around U = 7, a Mott gap opens in a continuous fash- 
ion (crossover). Since we cannot reliably study the low 
temperature behavior of the metallic phase within NCA, 
we will not investigate this transition in further detail. 
In the following, we will mostly focus on the insulating 
regime (U > 7). 

A good test of the implementation and its accuracy 
is the calculation of the total energy. The total energy, 
normalized by the number of sites in the transverse di- 
rection, has a local contribution 



N 



Epot — [U m d m + (ei oc , m — /i)r 



(28) 



where d m is the double occupancy and 7i m — Ti m ^ 

the occupation on layer m. In addition there is the intra- 

layer kinetic energy 



N 



^kin, intra ^ ^ ^ ^ ^h.m^h,m.o 



(29) 



m— 1 ha 

and the inter-layer kinetic energy 

JV-l 



£kin,intor = - ^ l ™ ( C k,m,* C k,m+l,<r) + h.C, (30) 



m— 1 ha 



where we have assumed vacuum boundary conditions. 



To evaluate Eq. (|30l) we note that t. 



-iti,G^^^ „j(M) and t m G k ,- 



;\ c k,m,a C k,m+l,o) — 
J k,m+l,m\ v i") ailKX <"m'~*to,m+l,m = ^h,m+l * Gh,m- 

The latter identity follows from a comparison of the 
Dyson equation §5§ with Eq. (fTTj|) . 

If an electric field is applied to the system, a current j 
will be induced (j is defined as the particle current, not 
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FIG. 3: Test of the energy calculation for a heterostructure 
composed of nine 1-d layers with t = 1, interaction U = 
15 on layers 1, 2, 3, 7, 8, 9 and U — 4 on layers 4, 5, 6, 
with "vacuum" boundary conditions. An electric field pulse 
is applied to all the layers. 



including the electric charge —1), 

3m ^ ^ (^fc^fc,m)^fc,m,CT; 
ha 

3m = ~i J^tnii^m.a^m+ha) ~ h.C. 



(31) 
(32) 



ka 



where j m is the intra- layer component, and j m is the 
current from layer m to m + 1. While the electric field 
is applied, the total energy will change like dE to t/dt = 
— '^2 im j m E m (for electrons with charge —1). Here we as- 
sume vacuum boundary conditions (and thus jg = jj^ — 
0), because otherwise energy can flow from the system 
into the leads. A good check of the numerics is thus to 
verify that E to t(t) — Ej(t) is time independent, where 

E j (t) = - />"[ ELi A (t) +EZZ{ 3 m (t)E m (t)] 
is the absorbed energy. After the pulse, the Hamiltonian 
of the system is time-independent, and the total energy 
should thus also become time-independent. In Fig. [3] we 
plot the time-evolution of the different energy contribu- 
tions for a nine-layer system consisting of three metallic 
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FIG. 4: Top panel: Time-dependent distribution of the dou- 
ble occupancy in a 39-layer system with U = 10, after a pulse 
excitation with fi « 12 applied to the middle layer. Lines 
show a fit to a Gaussian centered at the middle layer. The 
bottom left panel plots the widths extracted from such fits as 
a function of time for different values of the inter-layer hop- 
ping. The squared width grows linearly in t — t pu ise> where 
ipuisc ~ 1.7 is the time corresponding to the center of the 
pulse. Dashed lines show results with phonon bath (see text). 
The bottom right panel shows the dependence of the diffusion 
constant D on the inter-layer hopping. 



layers (U = 4) sandwiched between Mott insulating lay- 
ers (U = 15). The perturbation is an in-plane, few-cycle 
electric field pulse of frequency il « 12, which is applied 
to all nine layers. This strong pulse creates doublon-hole 
pairs and leads to a rapid increase in the potential energy. 
After the pulse, one observes a redistribution of potential 
energy into kinetic energy in such a way that the total 
energy is conserved. Also, the change in total energy 
is equal to the absorbed energy Ej, so that E to t — Ej 
remains zero within the numerical accuracy. That this 
result is a nontrivial check follows from the lower panels, 
which show the time-evolution of the double occupancy 



and intra-layer kinetic energy in all nine layers. These 
curves indicate that doublons and holes move from the 
insulating regions to the metallic region, where they re- 
combine, heat up the metal and lead to an increase in 
the intra-layer kinetic energy. 



B. Doublon diffusion 

As a first application we consider the spreading of 
photo-excited doublons and holes in a Mott insulator. 
The system consists of 39 layers and we employ the "re- 
peated" boundary condition to minimize boundary ef- 
fects. The doublons and holes are created in the central 
layer (m = 20) by the application of an in-plane electric 
field pulse with ft « 12, centered at t pu i sc — 1.7, which 
lasts up to t — 3. This set-up may not be realistic from an 
experimental point of view, but it allows us to study how 
artificially created carriers spread out inside a Mott insu- 
lating bulk. On the timescale of the present simulation, 
we can ignore the recombination of doublons and holes. 
This is consistent with corresponding DMFT calculations 
for a homogeneously excited bulk system, which indicate 
that the lifetime of these carriers depends exponentially 
on the interaction U in the Mott insulating regime J£ As 
shown in the top panel of Fig. 2] (results for U = 10), al- 
ready a short time after the pulse, the distribution of the 
photo-excited doublons (symbols) can be well fitted by a 
Gaussian (lines). For inter-layer hopping t x = 1, the 39- 
layer system allows us to track the motion of the doublons 
up to t » 20. Extracting the widths of the Gaussians 
and plotting them as a function of time (Fig. 01 lower 
left panel), we find that the square of the width grows 
proportional to t — t pu i so , indicating diffusive rather than 
ballistic motion. The doublon diffusion satisfies the ex- 
pected law d(m, t) - d cq (m, t) ~ exp(— (m — 20) 2 / (4-Dt)) 
(to = 20 is the central layer), with diffusion constant 
D m 1.03i x for t x = 1. As long as the doublon- holon re- 
combination is slow enough and the carriers are inserted 
with large kinetic energy, the diffusion of doublons and 
holes is not influenced much by the interaction strength. 
Within our numerical accuracy, we find the same diffu- 
sion constant for U = 7, 8, 9 and 10, even though U = 7 
is already close to the metal-insulator crossover. 

On the other hand, a smaller inter-layer hopping of 
course slows down the diffusion. The lower left panel of 
Fig-Hlplots the time-evolution of the squared width of the 
distribution, for t x ranging from 0.25 to 1.5 (U = 10). 
(Because of the rapid spreading of the charge carriers 
we cannot study much larger values of t x .) The diffu- 
sion constant D(i x ), which is extracted from linear fits 
to these curves, grows roughly quadratically with f x for 
small t x , while the dependence becomes almost linear for 
£ x > 0.5 (Fig.Hl lower right panel). 

In equilibrium, the diffusion constant is related to 
the conductivity a^c (with the charge e set to one) 
and the compressibility |^ via the Einstein relation (or 
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FIG. 5: Spreading of doublons in a 15-layer system with U = 
10 and pulse-excitation with Q « 12 on layer 1. The top panel 
shows the increase of the occupied part of the spectrum in 
layer 8, and the relaxation of the carriers towards the bottom 
of the upper band. The bottom panel shows the spectral 
functions on layers 1 and 8 (solid and dashed black lines) and 
the occupied part of the spectrum in layers 1, 5 and 8 a short 
time after the injection of carriers (rescaled in such a way that 
the maxima are approximately the same). The distribution 
of the fastest carriers remains almost unchanged from layer 1 
to 8. 



fluctuation-dissipation relation) 



a/j, 



Ode 



(33) 



Truly ballistic transport is thus expected for integrable 
one-dimensional systems (see Ref. |35| and references 
therein), which can have a perfect conductivity (i.e., a 
finite Drude weight <7dc ~ T>S(uj) for ui — > 0) even at tem- 
perature T > 0i2£ For the Hubbard model in higher di- 
mensions, the rather large width of the spectral function 
Ak a {uj) in the Mott insulator indicates that the scatter- 
ing time of a single particle excitation with momentum k 
is of the order of the inverse hopping, and hence its mean 
free path is not much larger than a few lattice spacings. 
Only in a Fermi liquid at T = would one expect infinite 



scattering times for electrons at the Fermi surface. 

To some extent, the behavior of D(t x ) shown in the 
lower right panel of Fig. U is qualitatively consistent with 
a quasi-equilibrium argument based on the Einstein rela- 
tion for large temperature T: Starting from the DMFT 
expression for the bulk conductivity^^ 



/oo 
. . _. - OO 



ka 



(a = j_, || ), the dc conductivity in the transverse direction 
and in the limit of high temperature is given by 

^dc«T^E«) 2 / (35) 

ka J -°° 

where v x = t x sm(k x ) is the band velocity perpendicular 
to the layers. The integral scales like 1/bandwidth, and 
the bandwidth is proportional to f" for i" 3> t x (almost 
independent layers), and proportional to t x for t x 3> i". 
Thus a x c ~ t x for t x > i 11 and a x c ~ \t x \ 2 /t" for i« > t x . 
Because |^ - n/AT for large T, the same behavior is 
found for the diffusion constant D[t x ). Physically, the 
behavior for small t x is consistent with a rate equation 
picture, where the transfer of a doublon from one layer to 
the next is given by Fermi's golden rule F ~ \t x \ 2 Af, with 
a matrix element oc t x , and a density of states J\f ~ 
that scales with the inverse bandwidth. 

Although the Einstein relation agrees with the ob- 
served behavior on a qualitative level, such a quasi- 
equilibrium theory cannot describe the spreading of dou- 
blons in detail. First of all, the initial perturbation of the 
system is strong, and it is neither clear on what timescale 
a local equilibrium description becomes possible, nor how 
well it would apply to a distribution that varies consid- 
erably over only a few lattice spacings. Since doublons 
and holes might cool down (lower their kinetic energy) 
while they spread in the bulk, equilibration could actu- 
ally lead to the formation of Fermi liquid quasi-particles 
and a corresponding reconstruction of the electronic den- 
sity of states, a process for which the time-scale is not 
known. Examples where nonequilibrium conditions have 
a strong influence on the spreading of particles have 
been studied recently, for a cloud of weakly-interacting 
ultra-cold atoms in an optical lattice (both fermions and 
bosons) i22r— For example, when the cloud expands into 
an empty lattice, it behaves diffusive in the dense core, 
but in the tails the density is too low to equilibrate, re- 
sulting in a ballistic expansion . 39 ' 40 

More detailed insight into the way in which doublons 
and holes spread into the bulk can be obtained from the 
time- and layer-dependent distribution function 



1 r°° 

A<(uj,t) = -Im dse^ s G<(t + . 
t Jo 



(36) 
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which reduces to the "photoemission spectrum" 
A < (u,t) — A(ui)f(uj) in equilibrium, and from the 
corresponding spectral function A m (w,t) (with G < 
replaced by —G R ). To study this quantity we switch to 
a smaller system, so that longer simulation times become 
possible and the integral in Eq. ([36]) does not strongly 
depend on to the upper cutoff. In the upper panel of 
Fig. [51 we plot the distribution function for a 15-layer 
system with U = 10, which is excited with a pulse with 
SI w 12 on the surface layer n = 1. (A "repeated" 
boundary condition is applied at layer 15.) On a given 
layer L, (L = 8 is plotted in the figure), the weight in the 
upper Hubbard band grows with time as more doublons 
arrive. At later times, the distribution is shifted to 
lower frequencies, indicating some kind of cooling of 
the particles as they move into the bulk. Still, the 
distribution is clearly non-thermal at all times, and its 
width remains comparable to the width of the Hubbard 
band. In such a highly excited system, one cannot 
expect the formation of quasi-particle states. Indeed, we 
only observe a slight broadening of the spectral function, 
rather than a formation of a quasi-particle band. 



Although the weight in the distribution function 
A^(uj,t) appears after an increasing-time delay as one 
moves further away from the surface, we find that the 
the distribution at the earliest times (i.e., right after it 
has achieved some measurable weight) has a similar shape 
on different layers (Fig.[5j lower panel). The distribution 
resembles the initial photo-doped distribution on layer 1, 
although the spectral function of the bulk layers is quite 
different from that of the surface layer, especially during 
the application of the pulse. This might be related to a 
coherent tunneling at early times. 



A detailed understanding of the various propagation 
effects at early and later times can be important to in- 
terpret the relaxation of photo-excited carrier distribu- 
tions in real experiments, which is governed by both dif- 
fusion and local relaxation phenomena. In real materials, 
doublons and holes can dissipate their energy to other 
degrees of freedom as they diffuse into the bulk, e.g., 
to phonons or spin excitations, which are not correctly 
accounted for in the DMFT formalism for the isolated 
Hubbard model. To study the consequences of this dissi- 
pation, we have simulated the diffusion in the presence of 
a local phonon bath with ojq = 1 and A = 1. In this case 
the doublons and holes spread more slowly, as shown for 
t x = 1 and t x = 0.5 by the dashed lines in Fig. [U A 
possible explanation is that the phonon cloud increases 
the effective mass of the carriers and hence reduces their 
diffusion coefficient. On the other hand, the curve for 
t x = 1.0 also reveals a slight negative curvature, which 
indicates that the cooling of the carriers influences the 
diffusion behavior in a nonlinear way. 
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FIG. 6: Top panel: Spectral functions for a heterostructure 
with five small-gap insulating layers (red, U = 10) on top 
of a bulk of large-gap insulator (blue, U = 20). The bands 
are shifted relative to each other in such a way that doublons 
can diffuse from the surface into the bulk, while holons can- 
not. The spectra are plotted with horizontal offsets of 0.3 
and arbitrary vertical offsets. Only four out of ten bulk lay- 
ers are shown. Bottom panels: Time evolution of the filling 
and intra-layer kinetic energy after excitation of the surface 
layer with a Q « 12 field-pulse on the surface layer. (Phonon 
coupling A = 1, phonon-frequency uio = 1, j3 = 10.) 



C. Surface excitation of a heterostructure and 
doping by diffusion 

An interesting application of the layer DMFT is to 
study the dynamics in heterostructures. Experimentally, 
such artificially designed systems may provide a way to 
confine the excitation to a well-defined region of the sam- 
ple (because, e.g., the pulse frequency can be tuned to 
the absorption band in certain layers), and induce con- 
trolled changes in the remaining layers. For illustration, 
we consider a heterostructure made of two different Mott 
insulators, and excite doublons and holes in the topmost 
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FIG. 7: The left panel shows the time-evolution of the den- 
sity, and the right panel the time-evolution of the double oc- 
cupancy in the set-up of Fig. [6] 



layer. As illustrated in the top panel of Fig. [13 the system 
consists of five Mott insulating surface layers (red spec- 
tral functions) on top of a Mott insulating bulk, whose 
gap is much larger than the gap of the surface layers (blue 
spectral functions) . The relative position of the Hubbard 
bands is chosen such that doublons can diffuse easily from 
the surface layers into the bulk, while the corresponding 
diffusion of holes into the bulk is prohibited. 

The diffusion of charge carriers leads to a time- 
dependent doping of the neighboring layers with elec- 
trons and holes, and the special setup of Fig. [5] allows 
to study the possible time-resolved emergence of a usual 
metallic state in the bulk layers, which are doped with 
electrons only. Explicitly, we simulated five surface layers 
with U — 10 on top of ten bulk layers with U = 20. We 
choose the "vacuum" boundary condition for the surface 
layer n = 1, and apply the in-plane electric field to this 
layer. To mimic dissipation to lattice and other degrees of 
freedom, which can accelerate the formation of a photo- 
doped state with low kinetic energy and less scattering, 
we couple the system to local phonon baths, as described 
in the methods section and in Ref. HU- The phonon bath 
parameters are ujq — 1 and A = 1. (The small structures 
visible in the spectral functions near the gap edges are a 
result of this phonon coupling.) 

The electron doping of the bulk- and net hole doping 
of the surface layers can be seen in the bottom left panel 
of Fig. |6l which plots the time-evolution of the density 
for the different layers. Note that even in the equilibrium 
system, a charge transfer occurs at the interface between 
surface and bulk layers, so that the first bulk layer is 0.2% 
electron doped, while the last surface layer is 0.2% hole 
doped. Figure [7] shows the time-evolution of the double 
occupancy and the density on a color scale (grey scale). 
Initially the double occupancy is slightly larger in the 



surface layers, due to the smaller value of U. We find that 
the interface between the two insulating regions does not 
slow down the diffusion of doublons into the bulk layers, 
while holes stay confined to the surface layers. There 
is even an accumulation of doublons on the bulk side of 
this interface, which is explained by a small downward 
shift of the Hubbard band. The net charge in the surface 
layers is reduced as time increases due to the holes which 
diffuse back from the interface. 

As a result of the dissipation, we expect the doublons 
and holes, which are created in the Hubbard bands of 
the L = 1 layer with a broad energy distribution, to cool 
down rapidly while they diffuse into the bulk. The latter 
effect should be evident as an accumulation of spectral 
weight in the distribution function (fM| at the lower edge 
of the upper Hubbard band (and symmetrically for the 
holes) . Figure|8]illustrates the time-evolution of the occu- 
pied spectral function in the upper Hubbard band, which 
roughly covers the energy range 1.5 < u> < 10. The few- 
cycle pulse with SI sa 12 creates doublons with a broad 
energy distribution centered at oj ~ 6 (in the middle of 
the upper band) . Such a broad spectrum is visible in the 
surface layer at t = 2 (the pulse lasts from about t = 0.4 
to t = 3). Very quickly (top right panel, t = 4), the dou- 
blons spread to the neighboring layers, and the cooling 
by the phonon bath leads to a shift of spectral weight 
to lower energies. Around t = 6, the diffusing doublons 
reach the bulk layers (n > 6). They keep diffusing into 
the bulk, which results in a pure electron doping of the 
bulk layers. Furthermore, by t — 10, the phonon bath 
has removed most of the excess kinetic energy so that 
the changes in the spectral function at later times are 
mainly due to changes in the carrier density. 

The decrease in the total in-plane kinetic energy in the 
different layers is also evident in the bottom right panel 
of Fig [6] This is consistent with a a metallization of 
the bulk layers as a result of the doping induced by the 
diffusion of doublons. The small quantitative change of 
the kinetic energy is explained by the small amount of 
doping and the high effective temperature of the doped 
system. Despite the strong coupling to the phonon bath 
with inverse temperature /3 = 10, the distribution func- 
tion A < (u),t) remains non-thermal within the accessible 
time-range, and it is much broader than expected for 
the effective temperature of the bath. In addition, no 
pronounced quasiparticle peak emerges in the spectral 
function on these timescales. As in the case of a photo- 
doped metallic state with electrons and holes^l it seems 
that the purely electron doped state obtained via dou- 
blon diffusion from the surface layer is not a good metal, 
and that the formation of a Fermi liquid state similar to 
an equilibrium chemically doped Mott insulator is a very 
slow process. 

Finally we note that in principle one should consider 
also the electrostatic energy associated with the (time- 
dependent) charge redistribution. This could be done 
for example by adding a layer-dependent Hartree poten- 
tial V m (t) to the chemical potential in the DMFT loop. 
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FIG. 8: Occupation function A < (uu,t) in the upper Hubbard band, for the same set-up as Fig. [15] The excitation is an in-plane 
field pulse with frequency Q m 12 on the surface layer. 



The formulas for this potential are given, for example, in 
Refs. Sill 

m— 1 k 

V m ({ni, . . . ,n m -l},t) — -a (m (t) — nbackground) i 

fe=l (=1 

(37) 

where ^background = 0.5 for half-filling and a is a constant 
proportional to the inverse dielectric constant. This po- 
tential would stop the spreading of charge into the bulk 
and confine the carriers to a region close to the inter- 
face. However, since the main purpose of the present 
work is to explain the nonequilibrium real-space DMFT 
method and to illustrate its versatility with several ex- 
amples, we will leave the calculation of realistic time- 
dependent charge profiles in heterostructures to a future 
publication. (The results shown here are representative 
of materials with a large dielectric constant.) 



D. Multi-layer structures under applied bias 

Transport through nanoscopic devices is another im- 
portant area of physics that involves both nonequilib- 
rium phenomena and strong correlations. The nonlin- 
ear current voltage characteristics of a two-terminal het- 
erostructures has been studied previously, using an inho- 
mogeneous steady-state DMFT approach. 43 The present 
formalism allows to study such systems in real time, and 
as a first application, we investigate the time-dependent 
build-up of current and charge distributions across the 
sample after the switch-on of a voltage-bias perpendic- 
ular to the layers. We consider a system consisting of 
L = 15 correlated layers in the Mott regime (U = 10). 
In these calculations we do not attach local heat baths, so 
that energy dissipation occurs only in the leads, and may 



not be relevant on the time-scales of our simulations. Ini- 
tially, the system is in equilibrium without applied bias, 
and at time t — 0, we switch on a bias v across the whole 
sample, assuming that the voltage drop is linear, i.e., the 
electric field is E ± = v/(N + 1). The top left panel of 
Fig. [9] shows the time-evolution of the current j 1 - flowing 
between the layers, for three different values of v. After 
some initial strong oscillations of the current j x , which 
are related to the build-up of a polarization perpendicu- 
lar to the layers, the currents into layer 1 (j'q ) and out 
of layer 15 (j^ 5 ) quickly settle to some w-dependent value 
which changes only slowly with time (bold lines). This is 
in contrast to the currents between layers in the interior 
of the sample, which show a slower time evolution and 
no relaxation into a quasi-steady state up to t — 30. The 
almost steady currents into and out of the leads exhibit 
a similar threshold behavior as was found in single-site 
DMFT calculations , 15 ! 16 i.e., an exponential increase at 
low bias of the form j 1 - oc v exp(— Vth/v). This is illus- 
trated in the top right panel of Fig. [5J which plots the 
ingoing and outgoing current at time t = 10 on a loga- 
rithmic scale. 

In the bottom left panel of Fig. Owe show current pro- 
files within the structure at different times, for v = 26. 
At short times, the current is largest near the leads and 
smallest in the center. Around t = 27, the current deficit 
in the center changes into a current surplus (see also up- 
per left panel), and we can expect some oscillations, until 
eventually an almost flat quasi-steady state distribution 
is established. This current profile implies a redistribu- 
tion of charge from the left side of the multi-layer struc- 
ture to the right side at short times. Indeed, a simi- 
lar plot of the density distribution (bottom right panel) 
shows a build-up of positive (negative) charge in the left 
(right) half of the structure which progresses from the 
boundaries. At t = 30 the excess charge peaks at layers 
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FIG. 9: Top left panel: Current between the different layers of a 15-layer system with U — 10 as a function of time after the 
switch-on of a voltage bias v. The black line shows the current which flows into the system, which at the resolution of the plot 
is indistinguishable from the current which flows out. Top right panel: Current measured at the leads at t — 10, divided by 
voltage and plotted as a function of voltage, showing an exponential increase at low bias. The slope gives the threshold voltage 
for the dielectric breakdown of the Mott insulator. Bottom panels: Current and charge distribution for different times. 



4 and 12, which is in the middle of the left and right re- 
gions. The distribution in the quasi-steady state might 
look similar. Again, one should in principle take the elec- 
trostatic potential associated with this charge redistribu- 
tion into account and compute the potential profile across 
the structure self-consistently. 

Similar time- and layer-dependent redistribution pro- 
cesses might be observable if they are triggered by a short 
pulse. To illustrate this, we finally discuss the current in- 
duced in Mott insulating structures under bias by an ap- 
plied intra-layer electric field pulse. We consider a 5-layer 
structure with U = 10. The voltage v = 2 across the in- 
sulating sample is small enough that after the build-up 
of a polarization, there is only a very small current flow- 
ing through the sample (Fig. [TU)) . Between t = 8 and 
t = 11.6 a field pulse with f2 « 12 is applied to the mid- 
dle layer (with polarization in the in-plane direction). At 



later times, the doublons and holes created by the pulse 
start to diffuse to the leads under the applied bias, which 
leads to a net negative current. The decay of this cur- 
rent is a direct measure for the mobilities. The intra-layer 
current during the pulse exhibits a peak in the opposite 
direction to the expected bias-induced current in the cen- 
tral region, which indicates that the polarization in the 
central layers is reduced in response to the perturbation. 



IV. CONCLUSIONS 

We have described and tested the nonequilibrium ex- 
tension of real-space DMFT, which allows to study lay- 
ered systems with strong electronic correlations. Like 
single-site DMFT (and in contrast to cluster-extensions 
of DMFT), the formalism is based on the assumption of 
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FIG. 10: Current induced in a Mott insulating 5-layer system 
(U — 10) under a voltage bias v — 2, by an in-plane field pulse 
acting on the middle layer. The pulse with f2 ps 12 is centered 
at fpuisc = 9.3. Initially, the system is in equilibrium, and the 
large current spikes at short times are due to the build-up of 
a polarization after the switch-on of the linear bias. 



a purely local self-energy. One thus only has to solve a 
collection of (coupled) single-site impurity problems in a 
self-consistent manner. For a layer geometry, in which all 
properties of the system depend on only one space direc- 
tion, the computational effort scales linearly with system 
size (up to the number of iterations, which may weakly 
depend on the system size), and the same is true for the 
storage requirement. We have discussed the details of our 
implementation based on self-consistent strong-coupling 
perturbation theory (NCA) as an impurity solver, but 
the formalism can equally be combined with a Monte 
Carlo, 44 or a perturbative weak-coupling solver 

As an application, we have simulated the diffusion of 
photo-excited doublons in a Mott insulator, both inside 
the bulk, and from the surface of a heterostructure into 
the bulk. The diffusion constant was found to depend 
mainly on the inter- and and intra-layer hopping, while 



it is almost independent of the interaction strength. A 
heterostructure set-up allows for a controlled doping of 
charge carriers of one type (e.g., doublons) into a Mott 
insulator, in contrast to photo-doping, where always both 
electrons and holes are inserted. In principle, this opens 
the possibility to study the formation of quasi-particles 
in a metallic system. For the current set-up, however, we 
find that the timescale for the build-up of such a state 
is rather long, such that the doped system behaves more 
like a bad metal on the numerically accessible timescales. 
A more thorough investigation of this important question 
will be deferred to a future study. 

The second type of application was the layer- and time- 
resolved calculation of the current through a correlated 
insulating slab, where we reproduced the threshold be- 
havior of the current-voltage characteristics known from 
previous nonequilibrium DMFT studies, and computed 
the evolution of the current- and density-profile after the 
switch-on of the voltage bias. We also studied a Mott 
insulating slab under bias (below the threshold for the 
dielectric breakdown) where the time-dependent redis- 
tribution of charge after a few-cycle laser pulse can be 
studied. 

In the future, one should include the effect of the elec- 
trostatic potential to obtain a more realistic description 
of the diffusion of electrons and holes in a heterostruc- 
ture. Also, the extension of our formalism to antiferro- 
magnetically ordered layers would be useful, because this 
would allow to exploit the cooling effect on the photo- 
doped carriers associated with demagnetization, 2 - How, 
and on which time-scale an almost thermal metallic state 
can be induced in a Mott insulator by diffusion of dou- 
blons from neighboring layers is an interesting topic for 
further studies. 
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